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We investigate the dynamics of a two-mode laser system by extending the two-mode Tavis- 
Cummings model with dissipative channels and incoherent pumping and by applying the mean-held 
approximation in the thermodynamic limit. To this end we analytically calculate up to four pos¬ 
sible non-equilibrium steady states (fixed points) and determine the corresponding complex phase 
diagram. Various possible phases are distinguished by the actual number of fixed points and their 
stability. In addition, we apply three time-delayed Pyragas feedback control schemes. Depending 
on the time delay and the strength of the control term this can lead to the stabilization of unstable 
fixed points or to the selection of a particular cavity mode that is macroscopically occupied. 
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I. INTRODUCTION 

Lasers build one of the key technologies in the current 
world as their rich dynamical behavior and high degree 
of control establish a solid basis for a wide range of appli¬ 
cations [1]. Especially, time-delayed feedback control [2] 
can effectively manipulate short and long time behavior 
of a laser system [3]. Typical examples are the control 
of laser bistability [4], chaos, and noise [5] as well as the 
manipulation of the laser emission [6, 7]. 

A common description of the controlled laser dynam¬ 
ics, particularly in the case of a quantum dot laser, is 
based on the senri-classical rate equations known as the 
Lang-Kobayaslri model [8]. It provides good agreement 
with the experiments if the photon output power is high 
enough [9]. However, there exists a more general mi¬ 
croscopic quantum treatment [10, 11], which describes 
successfully the photon statistics of laser light. It turned 
out that this microscopic laser theory also represents an 
essential ingredient for describing the Bose-Einstein con¬ 
densation of photons [12] which has been realized in dye 
filled microcavities in a seminal experiment in Bonn [13] 
and recently also in London [14]. Both lasing transition 
and Bose-Einstein condensation of light may appear in 
such systems under appropriate conditions, although the 
former reveals non-equilibrium physics, whereas the lat¬ 
ter represents an equilibrium phenomenon. For low cav¬ 
ity losses and above the external pumping threshold, the 
modes of the cavity become thermally populated accord¬ 
ing to a Bose-Einstein distribution with the macroscopi¬ 
cally occupied lowest mode [15]. However, for higher cav¬ 
ity losses the system behavior switches to be laser-like, 
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where one of the excited cavity modes becomes macro¬ 
scopically occupied and all thermal properties are lost 
[16]- 

Here we work out a two-mode laser model which al¬ 
lows to study under which conditions one of the two cav¬ 
ity modes becomes macroscopically occupied. To this 
end we extend the Tavis-Cummings model and consider 
N non-interacting two-level atoms in a two-mode opti¬ 
cal cavity with incoherent pumping and decay channels. 
Starting from a quantum master equation for the density 
operator we apply a mean-field approximation and deter¬ 
mine the equations of motion for the statistical averages 
of the respective system operators in the thermodynamic 
limit. We find an analytical solution for the steady states 
and obtain the resulting complex phase diagram. Under 
proper conditions, either the lower or the excited cav¬ 
ity mode can become macroscopically occupied. Hence, 
our model can be seen as a minimalistic precursor of the 
detailed model of photon condensation [12, 16]. In this 
sense, the former case could be referred as condensate-like 
and the latter case as laser-like state of light, although a 
direct analogy is not applicable due to the absence of the 
temperature scale in our simplified approach. The rich¬ 
ness of possible phases even within this reduced model 
indicates that the inclusion of realistic processes, like the 
thermalization via phonon dressing of the absorption and 
emission of the emitters, can potentially lead to an even 
larger variety of states. 

Additionally, we design different feedback control 
schemes to stabilize or to select one of the two radiat¬ 
ing modes. The two-mode laser, also known as two-color 
laser, with feedback was already studied both experimen¬ 
tally [17, 18] and theoretically [19]. However, these stud¬ 
ies within the Lang-Koboyashi model were focused on 
switching between the two modes using a non-Pyragas 
feedback type. In contrast to that we apply here the 
Pyragas type of feedback that was originally designed to 
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prevent chaos by stabilizing an unstable periodic orbit 
[20]. It is generally known as a powerful tool to change 
the stability of stationary states without modifying them. 
This is due to the fact that the feedback control term van¬ 
ishes in the stationary state since it is proportional to the 
difference of the system observable at two times t — r and 
t [ 21 , 22 ], 

The paper is structured as follows. In Sec. II we in¬ 
troduce the underlying model and apply a mean-field ap¬ 
proximation in the thermodynamic limit. In Sec. Ill we 
calculate the fixed points, investigate their stability, and 
discuss the resulting phase diagram. In Sec. IV we sug¬ 
gest several Pyragas feedback control schemes to stabilize 
the unstable mode or to select the mode of interest. Sec¬ 
tion V contains the summary of the obtained results with 
a short outlook. 


II. MODEL 

We consider TV non-interacting two-level atoms inside a 
two-mode cavity. The light-atom interaction is assumed 
to be of the Jaynes-Cummings type [23]. Thus, the total 
Hamiltonian of the system is 

2 2 

H = Widjdi + AJ Z + —j= J + + a\J ) (1) 

i— 1 * i —1 

and represents an extension of the Tavis-Cummings (TC) 
model [24, 25] from one to two modes. Here, we put 
h = 1 and (i € {1,2}) is a ladder algebra of the 

first/second cavity mode with frequency uqq, where we 
assume uq < u> 2 without loss of generality. The collective 
angular momentum operators are given by the sums J z = 
I Ef= i and = Ek= ! over all Pauli matrices of 
each two-level atom with energy level-splitting A. The 
population inversion of the atomic ensemble is directly 
related to J z , while its dipole moment can be expressed 
in terms of j . The coupling between the atoms and 
the optical mode assumes rotating wave approximation 
(RWA) and has the strength g/y/N that is taken to be 
the same for both modes. In spite of RWA, the TC model 
for large values of g has its own physical relevance since 
it can be experimentally realized in an ingenious setup 
using Raman transitions [26, 27]. 

To generate a lasing behavior and the interesting dy¬ 
namics we add decay channels and incoherent pumping 
to the system. We note in passing that two-mode Jaynes- 
Cummings models were studied in the past either with 
mode degeneracy [28, 29] or without dissipative effects 
[30] or without pumping of the atomic system but in pres¬ 
ence of additional driving of the cavity mode [31, 32]. 
Following Ref. [33], we couple our system to three dif¬ 
ferent baths. Both cavity fields are damped by coupling 
them to a zero temperature bath of harmonic modes with 
the characteristic decay rate re, while the atomic system 
radiates into the non-cavity modes with a rate yj_. Addi¬ 
tionally, the atomic system is incoherently pumped with 


a rate y-j-. Pumping can be formally described as cou¬ 
pling the atomic system to a bath of inverted harmonic 
oscillators [34]. All these effects are captured by the fol¬ 
lowing Markovian master equation of Lindblad type for 
the density operator p 

= ~ P\ ~ KL[a!]p - KL[d 2 ]p (2) 

N N 

-yE»-!EW' 

k =1 k =1 


with the Lindblad operator L[x\p = xA xp+pxJ x— 2xpxA. 
Pumping effectively occurs provided that y-j- > yj.. 

The dynamics of the statistical average (A) = Tr (Ap) 
of an arbitrary system operator A is described by 
d(A)/dt = Tr (Ap). To obtain a closed set of semi- 
classical equations, we perform the thermodynamic limit 
where the number TV of two-level atoms tends to infin¬ 
ity [35-39]. Therefore, we factorize the averages of an 
atomic operator A and a light operator L according to 
(alAj « (Al'j and rescale them with the atom num¬ 
ber TV, denoting the rescaled operator averages by cor¬ 
responding symbols without hat, i.e., J ± = (J±)/N, 
J z = ( J z )/N and a^l = (a^l)/y/N, where asterisk 
denotes complex conjugation. The resulting mean-field 
equations of the two-mode laser model are then 


(—re — zuq)ai — igJ ~, 

(3a) 

(—re + iu q)a{ + igJ + , 

(3b) 

(—re — iuq)a 2 — igJ~, 

(3c) 

(—re + ioj 2 )a 2 + igJ + , 

(3d) 

■ (-r_D - iA)J~ + 2ig(ai + a 2 )J z , 

(3e) 

■ (-r D + *A) J+ - 2 ig(al + a* 2 )J z , 

(3f) 


j z = T t (z 0 - J z ) + ig{a\ + a* 2 )J - ig(a\ + a 2 ) J + , 

(3g) 


where we have introduced the abbreviations Tt = ^Td = 

y; + Tr and z o = 2 ff t +^) ■ Note that J ~ = ( J+ )* and 

J, is a real quantity and by definition, one has —1/2 < 
zo < 1 / 2 . 

In the one-mode limit, the corresponding equa¬ 
tions similar to Eqs. (3) represent a common ex¬ 
ample of a laser model. For the critical value of 


‘2 j 2 

< 7 C = {1 + optical mode becomes 

macroscopically occupied, i.e., a phase transition occurs 
from a non-lasing to a lasing state. In the limit of van¬ 
ishing pumping and losses, i.e. Fp—>0 ,k—> 0, Eqs. 
(3) describe the quantum phase transition in the Dicke 
model with RWA from a normal to a superradiant phase 
[37, 40-43]. Thus, the presence of the two modes and 
the pumping term allows the generation of a much more 
complicated dynamics, as either of the two modes can 
be macroscopically occupied. Moreover, we can influence 
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the dynamical evolution of the system by applying differ¬ 
ent Pyragas time delay schemes, which allows to stabilize 
or destabilize the modes and to select the transition type. 


III. DYNAMICS WITHOUT FEEDBACK 

Equations (3) describe the dynamical evolution of the 
two-mode system depending on decay rates and pumping 
strength. Steady state of these equations can be either a 
stable fixed point or an oscillating state , i.e. a limit cycle. 
In the following we provide an analytical description of 
the possible steady states. 


A. Steady states 

The system (3) has a trivial fixed point a® = «2 = 
K)° = (a$)° = 0, (J+)° = (J-)° = 0, and J° = z 0 , 
where no cavity mode is occupied and the atomic en¬ 
semble has a stationary population inversion with zero 
dipole moment. Due to the 17(1) symmetry of the Eqs. 
(3), there also exist non-trivial solutions that can oscil¬ 
late in time with some characteristic frequency, so that 
the observables, like the mode occupation a*ai, reach a 
fixed value. To find such steady state solutions, we have 
to determine the frame where also a\ 2 , and J ± reach a 
fixed value. Therefore, we switch into a frame rotating 
with frequency w, which has to be determined, i.e. we 
put at -a aie~ iu,t , a* -A a*e iuit , J± -A J ± e ± ^ t . Note 
that this transformation shifts the natural frequencies of 
both the cavity modes and the atoms by w, i.e., 

w-i —>■ Wi — to = w,; jS , A -a A — u) = A s , (4) 

but does not change the observables like a}ai. Setting 
d ^*2 in the transformed Eqs. (3a-d) to zero, we can ex¬ 
press these cavity quantities in terms of J ± . Next, setting 
./^ to zero in the transformed Eqs. (3e,f) with the cavity 
quantities being eliminated, we find the requirement 

0 = J ± | ± 2g 2 J z [ 2k + i(wi iS + W 2 jS )] 

T (T^i T tA s )(/v T iuj\ s )(k T iuj 2iS ) j'■ 

! 

For J ± ^ 0 the previous equation determines the 
of the stationary atomic inversion 

j 0 _ (Tj - ?A s )(k - ZUTs)(K - iV2,s) 

2g 2 (2n - zwi jS — zw 2 , s ) 

However, since J z has to be real on physical grounds, its 
imaginary part has to be zero. This condition enforces 
the characteristic frequency co to solve the equation 


(5) 

value 

( 6 ) 


Td(wi, S + w 2,s) (k 2 + Wi jS W 2 >s ) 

+kA s (2k 2 + uj 2 s + u >2 S ) = 0. (7) 


Note that, due to Eq. (4), Eq. (7) is a cubic equation in 
uj and has up to 3 real solutions. For each real solution 
w, the real part of the expression for J z in ( 6 ) gives the 
steady state expectation value 


k(T 2 d + A}) (2k 2 + u;f s + w 2 )S ) 
2g 2 T D [4k 2 + (wi, s -t-w 2jS ) 2 ] 


The remaining transformed equation (3g) can be solved 
for J + J~ in the steady state, yielding 


(J+J-)° 


r T (z 0 - J°) (k 2 + uls) (k 2 + tufj 
2g 2 K (2k 2 + ul s + 


(9) 


Since J + J~ has to be positive, the obtained steady state 
values are physical iff J z < zq. If that is the case, the 
previous equation fixes up to the phase factor. There¬ 
fore, we may choose (J + )° = (J _ )° = \J (as a 
steady state expectation. Finally, the corresponding ex¬ 
pressions for o° and (a *) 0 (i G { 1 , 2 }) in terms of (J *) 0 
follow from their transformed equations 


o _ ig(J )° _ ig(J + )° 

'i | • i ) 

K + ZUAj ;S K — ZWj iS 


( 10 ) 


With this we have found a complete set of steady state 
solutions for our two-mode model. Each physical solution 
for a characteristic frequency w corresponds to a different 
non-trivial fixed point. Thus, together with the trivial 
fixed point, the laser model possesses up to 4 different 
steady state configurations, whose stability properties we 
are going to study in more detail in the next subsection. 


B. Stability of steady states 

First, we investigate the stability of the fixed points. 
This is checked as usual by linearizing the mean-field 
equations (3) in the rotated frame around the fixed 
point and by determining the eigenvalues of the lin¬ 
earized system. An eigenvalue with positive (negative) 
real part would support the solution divergence (conver¬ 
gence) from (to) the fixed point, which is then unstable 
(stable). If not mentioned otherwise, we choose the fol¬ 
lowing parameter values: oj\ = 2A, w 2 = 4A, 74 . = 0.1A, 
7 t = 0.2A. 

Fig. 1 shows the main result in form of a complex 
phase diagram in the g-K plane for two different pumping 
rates 74 - = 0.2A,0.5A, encoding the total number and 
the number of stable fixed points. We see that, if the 
atom-field coupling is too small, only one trivial solution 
exists which corresponds to region (a). By overcoming 
some critical value for g, at least one non-trivial solu¬ 
tion appears, thus the wi and W 2 modes become rnacro- 
scopically occupied. For smaller K-rates, we see a rich 
structure in the phase diagram. One can have different 
combinations of possible and stable fixed points, which 
are represented by a combination of color and dashing 
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FIG. 1. (Color online) The phase diagram shows the total 
number of fixed points and the number of stable fixed points 
in the g-K plane. For small k, there exist up to 4 physical fixed 
points, 2 of which are stable. In the region (c) all fixed points 
are unstable. Table I sums up the main properties of the re¬ 
gions (a)-(g). The green color gradient encodes the mode pop¬ 
ulation ratio ni/ 712 , where n; = a*at. The lower part shows 
the effect of increased pumping. Parameters: uj\ = 2A, 0 J 2 = 
4 A, 74 . = 0.1A, 7 -f = 0.2A (upper), 7 ^ = 0.5A (lower). 


in Fig. 1. For example, the region (d) has two non¬ 
trivial physical solutions, but only one is stable. Table 
I provides the corresponding overview. For larger k and 
^-values, the phase diagram contains region (c) without 
any stable fixed points. Here the system observables, like 
the mode occupation, oscillate with fixed frequency and 
amplitude, thus a limit cycle represents the only stable 
solution in this area. Note, that we have found no sta¬ 
ble limit cycles except in region (c). The coloring in the 
(b)-region shows the ratio ni/ri 2 of occupation of both 
modes, where rii = a*a,i. We observe that the occupa¬ 
tion ratio and thus the dominating mode changes with 
the dissipation rate k and the coupling strength g. Note 
that in the regions (e) and (g), where we have two sta¬ 
ble fixed points, both ratios n\/ri 2 ^ 1 for fixed n and g 
values exist. Especially in this region one of the modes 
is much more occupied and vice versa, thus the emitted 


Area 

(a) 

(b) 

(c) 

(d) 

(e) 

(f) 

(g) 

#(FP) 

1 

2 

2 

3 

3 

4 

4 

#(SFP) 

1 

1 

0 

1 

2 

1 

2 


TABLE I. Overview of the total number of fixed points #(FP) 
and the number of stable fixed points #(SFP) within different 
regions of the phase diagram in Fig. 1. 


radiation comes here mainly from one mode. 

The lower part of Fig. 1 shows the effect of increased 
pumping. We see that the region with more than two 
fixed points (d)-(g) becomes larger, while the limit cycle 
region (c) is shifted to higher k values. 

Fig. 2 shows the occupation of both modes as a func¬ 
tion of coupling strength g for a fixed value of n = 0.01 A, 
along the horizontal grey arrow in the phase diagram of 
Fig. 1. We plot all possible stationary solutions including 
the unstable ones. The unstable fixed points are dashed, 
the occupations, which belong to the same fixed point, 
have the same color and the same thickness. The curves 
of the second mode are additionally marked with crosses. 
We see different types of bifurcations while increasing g. 
First, at g = 0.3A a pitchfork bifurcation occurs, where 
the trivial solution becomes unstable and a new stable 
solution occurs. Afterwards, an additional bifurcation 
takes place at g = A, where an unstable solution splits 
up from the trivial one and becomes stable at g = 1.5A. 
Later, at g = 3.2A, a third bifurcation with an unstable 
solution splits up. For the used parameter values Eq. (7) 
has three real roots, nevertheless at least one of the ob¬ 
servables in Eqs. (8)-(9) is unphysical, for instance a neg¬ 
ative mode population rii or an imaginary J + J~ value. 
Thus we have only two non-trivial solutions for g > 1.5A. 
The two solutions allow the lower or the upper mode to 
have a high occupation, respectively. Note that the so¬ 
lution depends crucially on the chosen initial condition. 
Fig. 3 shows an example of this behavior where we vary 
the initial state of the cavity modes ni(0), n 2 (0) for a 
given initial state of the atomic system. In the light blue 
area (diagonal lines) the system converges to the fixed 
point FP 1, in the dark blue area (vertical lines) to the 
fixed point FP 2 from Fig. 2. 

In the next section we present different Pyragas feed¬ 
back schemes. They allow to switch between a macro¬ 
scopic occupation of the two cavity modes irrespective 
of the chosen initial condition and also to change fur¬ 
ther dynamical properties like the fixed point attraction 
region of the considered model. 


IV. DYNAMICS WITH FEEDBACK 

We now demonstrate the impact of time-delayed feed¬ 
back control on the system. As a feedback signal we 
always use one of the system properties and restrict our¬ 
selves only to Pyragas feedback type [20]. Therefore, we 
insert into the mean-field equations Eq. (3) an additional 
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FIG. 2. (Color online) All stationary solutions of the mean- 
field Eq. (3) for the occupation of both modes (m,ri 2 ) are 
plotted as a function of g for fixed K-value along the horizon¬ 
tal dashed arrow in Fig. 1 (upper). The unstable solutions 
are dashed, the solution set is marked by the same color and 
the same thickness. The trivial solution with zero-mode oc¬ 
cupation is always present but unstable beyond a critical g. 
Note, that all occupations in the plot are shifted by 10 ~ 2 due 
to the log-scaling. Parameters: k = 0.01 A, on = 2 A,W 2 = 
4A,7j_ = 0.1A,7-|- = 0.2A. 


Dominant 

attraction: 



A. Stabilization of fixed points 

The phase diagram in Fig. 1 has regions with non¬ 
trivial unstable steady states, which do not attract the 
solution. If no stable point exists, the solution oscil¬ 
lates periodically. This occurs only in the region (c), 
see grey dotted curve in Fig. 4 (left) obtained using the 
parameters k = 0.5A, g = 5A. To stabilize the unstable 
non-trivial fixed point we suggest the following feedback 
scheme of Pyragas type [20] 

j z ^j z - A [J z (t - t) - J z {t )\, (11) 

thus we modify the population inversion by a difference of 
the J z spin component at two different times t — r and f, 
where r denotes the time delay parameter. Additionally, 
this difference is scaled by A. The feedback term in Eq. 
(11) can be realized, for instance, by extra pumping of the 
atomic system or by opening additional decay channels, 
depending on the value of the feedback signal A [J z (t — 
T)-Ut)]. 

The solid lines in Fig. 4 (left) show feedback actions 
for a point in the region (c). We see, that for t 1/A 
the mode occupations become constant, thus the fixed 
point is stabilized and the feedback signal vanishes. In 
contrast, without feedback the oscillations with finite am¬ 
plitude are always present (grey dotted line). The right 
part of Fig. 4 shows the control diagram [44] in the r-A 
plane. The color encodes the largest real part of all ex¬ 
isting eigenvalues, obtained from the linearized equation 
of motion [21] (see Appendix Al). The fixed point is 
stable if this value is negative, which is the case in the 
blue area (Fig. 4, right). For the boundaries (green dots 
in Fig. 4, right) an analytical expression can be derived, 
see Appendix A1. 


FIG. 3. (Color online) Attraction region of two stable fixed 
points from Fig. 2 depending on the initial population of the 
cavity modes ni(0) and ri 2 ( 0 ). Used parameters: J + (0) = 
J-(0) = 0.185, J z { 0) = 0.076, g = 2A, k = 0.01A, wi = 2A, 
0 J 2 = 4A, 7 j. = 0.1A, 74 - = 0.2A. 

control term, which is conditioned on the difference of a 
system property at two different times t — r and t, where 
r represents a time delay between the signal determina¬ 
tion and the feedback into the system. Due to the rich 
phase diagram even without feedback in Fig. 1, it seems 
impossible to engineer one feedback scheme which works 
in every part of the phase diagram. Hence, we have to 
find for each part of the phase diagram a scheme which 
produces the desired results like mode selection or sta¬ 
bilization. However, the chosen feedback may not work 
in other parts of the phase diagram or will have other 
influences onto the system dynamics. In the following, 
we present three feedback schemes for different purposes 
and parts of the phase diagram, give a possible imple¬ 
mentation picture for each scheme and demonstrate ex- 
emplarily their influence onto the system evolution. 


B. Selection of the dominantly occupied mode 

We now focus onto the region (e), which features two 
stable non-trivial fixed points. The main interest in this 
region is the occupation of the respective cavity modes. 
In each of both solutions one mode has a high occupa¬ 
tion, whereas the other one has a low occupation, see 
Fig. 2. In that way, the light leaking out from a cavity 
is generated by mostly one of the two modes. Without 
feedback the dominating mode is selected by the initial 
condition, see Fig. 3, which is usually hard to control. In¬ 
terestingly, we found a feedback scheme, which allows to 
select the mode of interest, i.e. to select the frequency of 
the radiated light, which was also achieved for a quantum 
dot laser in the Ref. [19] with a non-Pyragas feedback 
type. We argue that our feedback type can switch the 
system behaviour between a macroscopic occupation of 
the higher or the lower cavity mode. 

To select the lower mode uq we modify its frequency 
in Eqs. (3) as 

u>i -a wi + A [n 2 (t - t) - n 2 {t )], 


( 12 ) 
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FIG. 4. (Color online) (left) Pyragas feedback control of J z 
( 11 ) stabilizes the non-trivial fixed point in region (c) of phase 
diagram Fig. 1. Without feedback the stationary solution is 
a limit cycle (gray dotted curves). With feedback the solu¬ 
tion converges to a fixed point (solid curves). Parameters: 
r = A -1 , A = 0.4A. (right) Control diagram in r-A-plane. 
Vertical scale bar gives the largest real part of the eigenval¬ 
ues of the linearized equations. In blue region fixed point 
becomes stable. Green dots show the boundaries from an 
analytical expression, see Eq. (A7). Parameters: re = 0.5A, 
g = 5A, wi = 2A, ui 2 = 4A, 74 , = 0.1A, 74 . = 0.2A. 


where ri 2 = a^a 2 represents the occupation of the second 
mode. This feedback type is also measurement based as 
the mean photon flux is proportional to the mean oc¬ 
cupation of the photonic modes [45, 46]. Thus, the fre¬ 
quency of the first mode has to be changed according to 
the difference of mean photon fluxes of the second mode 
at times t — r and t. 

However, the previous (or similar) feedback scheme 
does not work well for selecting the upper mode u> 2 . For 
that purpose we modify the feedback scheme according 
to [47] 

di —> a\ — A[ai(f — r) — ai(f)], (13) 

which is now a coherent type of feedback, as one can in¬ 
terpret it as a direct control without measurement [47]. 
One possible realization is the back coupling of emitted 
photons by a mirror, where the mirror distance fixes the 
time delay r [48]. This scheme works for a properly cho¬ 
sen T-parameter [21] as, for instance, r = 27r/w (or mul¬ 
tiples of it), where u> denotes the characteristic frequency 
of the rotated frame determined by Eq. (7). This choice 
guarantees that the feedback term in Eq. (13) vanishes 
for t 1/A. 

The action of both feedback types is shown in Fig. 5 
for the system parameters re = 0.005A,g = 2A and the 
feedback parameters A = 0.01A,t = A -1 (upper) or 
A = A,t = 2tt/oj (lower), where w denotes the rotating 
frame frequency determined from Eq. (7). Solid marked 
curves show the cavity mode occupations with feedback, 
dashed curves without feedback. Both feedback schemes 
destabilize only one fixed point in the region (e) of Fig. 
1, thus the system converges to the other one. In the top 
figure we see the action of feedback Eq. (12). Without the 
feedback, the excited mode ui 2 has a dominant population 
(dashed violet line), whereas with control its occupation 



FIG. 5. (Color online) Usage of feedback schemes in region (e) 
of Fig. 1 for driving the system toward a macroscopic occupa¬ 
tion of the lower (top) or higher cavity mode (bottom), (top) 
Feedback scheme Eq. (12) selects highly populated ground 
mode (red line with markers), whereas (bottom) control type 
Eq. (13) selects highly populated excited mode (violet line 
with markers). The inset (bottom) shows the zoom for small 
photon numbers. Without feedback the other modes have a 
macroscopic population (dashed violet and red lines in both 
figures). Parameters: A = 0.01A (top), A = A (bottom), 
re = 0.005A, g = 2A, wi = 2A, ui 2 = 4A, 74 . = 0.1A, 
7 t = 0.2A. 


becomes low (violet line with markers) and instead the 
ground mode u> 1 (red line with markers) is macroscop- 
ically occupied. The bottom figure shows the opposite 
behavior. Instead of the lower mode (red, dashed), the 
higher mode is macroscopically occupied (violet line with 
markers). Note, that both stable steady states exist with¬ 
out feedback in the region (e) of Fig. 1. However, their 
attraction regions depend on the initial condition, as is 
shown without feedback in Fig. 3. We emphasize, that 
with feedback the selection of modes works independently 
of the chosen initial condition for the tested parameter 
values. 

Fig. 6 shows the control diagram in the r-A space with 
re = 0.005A, g = 2A for the feedback type Eq. (12) ob- 
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tained from a linear stability analysis. We see that there 
are parameter regions where only one of the fixed points 
becomes unstable and also where both fixed points be¬ 
come unstable. In the blue-dotted area the fixed point 
with 7i2 ^ becomes unstable, whereas in the green- 
dashed region another fixed point with m n ,2 is desta¬ 
bilized. The boundaries are calculated analytically (see 
Appendix A 2). In order to reach the fixed point with a 
macroscopic occupation of the lower cavity mode we have 
to choose the parameters in the region having only blue 
dots. Fixing the feedback parameter in the region hav¬ 
ing only green dashes (arrow in the diagram) should se¬ 
lect the fixed point with a macroscopic population of the 
higher cavity mode. However, there are some exceptions. 
The fixed point with ri 2 7ii attracts the solution if the 
initial condition is rather close to it, otherwise the solu¬ 
tion converges to a limit cycle which appears in this case 
in presence of Pyragas control [22]. Limit cycle solutions 
are also present in the parameter area where both fixed 
points become unstable due to the time-delayed feedback 
control. 



FIG. 6. (Color online) Stability diagram for Pyragas feedback 
type Eq. (12). In the dashed (dotted) region the first (second) 
fixed point (FP), related to a macroscopic population of the 
lower (higher) cavity mode as in Fig. 2, becomes unstable. 
Parameters: k = 0.005A, g = 2A, uj\ = 2A, u )2 = 4A, 
74. = 0.1 A, 7 t = 0.2A. 


V. SUMMARY AND OUTLOOK 

In this paper we have investigated the mean-field dy¬ 
namics of a two-mode laser model based on an ex¬ 
tended Tavis-Cummings model in the thermodynamic 
limit without and with time-delayed feedback. The corre¬ 
sponding mean-field equations can be solved analytically 
in the steady state. Even without feedback control this 
model exhibits a complex phase diagram with multiple 
stable fixed points. Our Pyragas feedback schemes allow 
to drive the system to different phases, by selecting or 


stabilizing one preferred stationary solution. 

We studied also other feedback schemes of the Pyra¬ 
gas type, but they led to similar results as already shown. 
However, especially in phases with a combination of un¬ 
stable and stable non-trivial fixed points it is difficult 
to design a feedback scheme which stabilizes or selects 
one stable configuration for a wide range of initial con¬ 
ditions. The reason for this is that the Pyragas control 
type affects the stability of all fixed points. For example, 
the stabilization succeeds only close to the correspond¬ 
ing fixed point in the sense of the linear stability anal¬ 
ysis. Farther away from the fixed point, we have often 
observed the appearance of limit cycles with large attrac¬ 
tion regions or even chaotic solutions, which is a known 
feature in laser systems with feedback [49] and also oc¬ 
curs for other non-linear dynamical systems with time 
delay [50-53]. 

Since our calculations were done at a semi-classical 
level by restricting ourselves to first-order cumulants, we 
expect that the results should hold in the thermodynamic 
limit, where the number N of two-level atoms tends to in¬ 
finity. On the one hand, the fluctuations scale like 1/y/N 
with the number of atoms N [34]. On the other hand, 
the laser dynamics or a condensation is usually studied 
at this level. Furthermore, the semi-classical regime of 
the quantum-optical models like Dicke [54] or Lipkin- 
Meshkov-Glick [55] predicts correctly their main proper¬ 
ties, like observable averages or occurrence of a quantum 
phase transition [37, 39, 56]. However, going beyond the 
factorization assumption could be performed by includ¬ 
ing higher-order cumulants, e.g. by using the Gaussian 
approximation, which involves first- and second-order cu¬ 
mulants [36, 38, 57]. 

It would be certainly interesting to analyse the impact 
of control on the quantum fluctuations. This could be 
investigated with other approaches to feedback [58, 59], 
which usually requires a high numerical effort. In this 
respect a promising feedback scheme was introduced in 
Refs. [60, 61], which allows to control the entanglement 
and light bunching by structured environment and con¬ 
verges to a Pyragas control type in the one excitation 
limit. However, the general quantum version of Pyra¬ 
gas control type remains an unsolved question. A new, 
conceptually significant approach has been recently intro¬ 
duced in Ref. [62], although it appears to be numerically 
demanding. 

Finally, we note that it would be worthwhile to ex¬ 
tend our two-nrode laser model with the thermalization 
mechanism along the lines of Refs. [12, 16]. This would 
yield a minimal model to study the transition between 
a condensate- and a laser-like state which originate from 
a macroscopic occupation of the lower and higher cav¬ 
ity mode, respectively. Adding Pyragas feedback control 
terms as suggested here should, thus, allow to switch the 
system behaviour between condensate- and laser-like. 
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Appendix A 

In the following we show how to determine the bound¬ 
ary condition in the stability diagrams Fig. 4 (right) and 


of Fig. 6 in the presence of time-delayed Pyragas feedback 
control terms Eq. (11) and Eq. (12), respectively. 

1. Stabilization of fixed points 

Linearizing the equation of motion (3) together with 
the feedback condition Eq. (11) we obtain the equation 

<5v(t) = A Sv(t) + B Sv(t — t), (Al) 

where v = (ai, a*, a 2 , a 2 , J + , J - , J z ), <5v gives a devia¬ 
tion from the fixed point v° determined via the proce¬ 
dure given in Section III A and we have introduced the 
matrices 


—ioj ! jS — K 

0 

0 

0 

0 

-ig 

0 

0 

iiUl'S — K 

0 

0 

ig 

0 

0 

0 

0 

—i 0 J 2 ,a — K 

0 

0 

-ig 

0 

0 

0 

0 

— n 

ig 

0 

0 

0 

-2 igJ° z 

0 

-2 igj° z 

iA s — To 

0 


2 igJ° 

0 

2 igj° 

0 

0 

— iAs — Td 

2 ig{a°i + a 2 ) 

-ig(J + )° 



ig(J~)° 

-ig{a\ + a°) 

ig((a*)° + K) 0 ) 

-r T 


B = A • (0,0,0,0,0,0,1) T ( 8 ) (0,0,0,0,0,0,1). 


The stability condition is then [21] 


where 


0 = det [(A-B)-B- e~ Ar - Al]. (A2) 

The fixed point is stable if all possible solutions for A 
have negative real part. From Eq. (A2) the equation 
for phase boundaries can be obtained as follows. At the 
phase boundaries, A has vanishing real part. Thus, re¬ 
placing A —> ifl (12 £ R) in Eq. (A2) and calculating the 
determinant, we obtain 


6 7 


'” >>.,,107 . y ri/) o. 

(A3) 

0 

II 

■Oj 

0 

II 

■05 


I 1, j even, 

1 i, j odd, 

(A4) 


Cl = B 0 + B 2 n 2 + B 4 n 4 + B e n 6 , (A6) 

c 2 = a 0 + a 2 q 2 + a 4 d 4 + A 6 n 6 , 
c 3 = a 1 d + a 3 o 3 + a 5 o 5 , 
c 4 = Bid + B 3 n 3 + B 5 n 5 + B 7 n 7 . 


Squaring and summing the Eqs. (A5), we can eliminate 
the r dependence and obtain a 14th order polynomial 
equation in D. This provides up to 14 solutions for D, 
but only two of them turn out to be real. Next, we sum 
both of the Eqs. (A5) together in a suitable way in order 
to eliminate the sin term. The resulting equation can 
then be solved for r as 


where A i: Bi , i £ (1, 2,... .7} are real coefficients which 
depend on the system parameters both explicitly and im¬ 
plicitly via the fixed point solution, and on the feedback 
strength A. However, the corresponding expressions are 
too long for showing them here. 

Splitting the equation in real and imaginary part, we 
obtain the following two equations 

0 = Gi + C 2 cos(flr) + C3 sin(f2r), (A5) 
0 = C4 + C3 cos(Dt) — C 2 sin(flr), 


1 ( C3C4 + C 4 C 2 \ ,277 ^ ^ 

7 = a (— Cj + Cj j + 77-’’ * £ z - <A7) 

This yields the boundaries in Fig. 4 (right), which per¬ 
fectly agree with the corresponding numerical calcula¬ 
tions. Two valid solutions for O build the IJ-shaped 
structure in the diagram, whereas 2 is responsible for 
its periodic structure. 
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2 . Selecting the fixed point Eq. (A3) 

4 7 

0 = e~ inT £ Cj A. ,-SV + ]T c,!},<!'. (A8) 

j=0 j=0 

The procedure is similar to Appendix A 1, but the feed¬ 
back condition is given now by the Eq. (12). The matrix 
B is then redefined as 

As the parameters A :j . Bj are real, Eq. (A8) can be split 
in real and imaginary parts, which yields 


1 , j even, 
i, j odd. 


(A9) 


B = -iX 


(0 0 
0 0 
0 0 
0 0 
0 0 
0 0 
v o 0 


(a$)V° 

-a 2 °K)° 

0 

0 

0 

0 

0 


{a* 2 )°a i° 0 0 0\ 
-a 2 0 (Aj;) 0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 Oy 


0 = C\ + C 2 cos(rir) + C 3 sin(fir), (A10) 

0 = C 4 + C 3 cos (fir) — C 2 sin(flr), 

where 

Ci = B 0 + B 2 ^ 2 + B 4 n 4 + B 6 n 6 , (All) 
c 2 = A 0 + A 2 n 2 + A 4 n 4 , 

C 3 — .111 1 T- 

c 4 = b x vl + B 3 n 3 + B 5 n 5 + B 7 n 7 . 


From the upper equations one can then eliminate the 
r dependence to determine possible 41 values. With this 
r can be calculated as in Eq. (A7), but 6); is then re- 
The further procedure is the same. First we calculate the placed by C). The resulting (42, t) combinations are the 

determinant Eq. (A2) and write it in a similar form of boundaries in Fig. 6. 
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